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Abstract 



This paper studies state quantization schemes for feedback stabilization of control systems with 
limited information. The focus is on designing the least destabilizing quantizer subject to a given 
information constraint. We explore several ways of measuring the destabilizing effect of a quantizer 
on the closed-loop system, including (but not limited to) the worst-case quantization error. In each 
case, we show how quantizer design can be naturally reduced to a version of the so-called multiccnter 
problem from locational optimization. Algorithms for solving such problems are discussed. In partic- 
' ular, an iterative solver is developed for a novel weighted multicenter problem which most accurately 

D"""- . represents the least destabilizing quantizer design. 

O 

m 



1 Introduction 



In this paper we study control systems whose state variables are quantized. We think of a quantizer as 
a device that converts a real-valued signal into a piecewise constant one taking a finite set of values. 
The recent papers g g, H discuss various situations where this type of quantization arises and provide 
references to the literature. Mathematically, a quantizer can be described by a piecewise constant 
. £h ' function q : T> C W 1 — > Q, where Q is a finite subset of R n with a fixed number of elements N. Here 
^ . n is the state dimension of a given system. We denote the elements of Q by qi, . . . , and refer to 
^ | them as quantization points. The sets Wi := {x G V : q(x) = i G {1, ... , iV} associated with fixed 
values of the quantizer form a partition^ of the domain T> and are called quantization regions. Quantized 
values of the state represent a limited information flow from the system to a feedback controller: the 
state is not completely known to the controller, but it is only known which one of a fixed number of 
quantization regions contains the current state at each instant of time. 

In the literature it is usually assumed that quantization regions are fixed in advance and have specific 
shapes, most often rectilinear. Here we are interested in the situation where the number of values N 
of the quantizer is a given information constraint in the control problem, but the control designer has 
flexibility in choosing a specific configuration of quantization regions (whose shapes can in principle be 
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1 A collection {Wi, . . . , Wn} of subsets of I? is a partition of T> if the intersection between the relative interior of any 
two Wi is empty and the union of all Wi equals T>. 
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arbitrary) and quantization points. While there has been some research on systems with quantization 



regions of arbitrary shapes [14, 12 1 and on the relationship between the choice of quantization regions 
and the behavior of the closed-loop system (?], the general problem of determining the "best" 
quantizer for a particular control task such as feedback stabilization remains largely open. 

A feedback law which globally asymptotically stabilizes a given system in the absence of quantization 
will in general fail to provide global asymptotic stability of the closed-loop system that arises in the 
presence of state quantization. There are two phenomena which account for changes in the system's 
behavior caused by quantization. The first one is saturation: if the quantized signal is outside the 
range of the quantizer, then the quantization error is large, and the control law designed for the ideal 
case of no quantization leads to instability. The second one is deterioration of performance near the 
equilibrium: as the difference between the current and the desired values of the state becomes small, 
higher precision is required, and so in the presence of quantization errors asymptotic convergence is 
impossible. These phenomena manifest themselves in the existence of two nested invariant regions such 
that all trajectories of the quantized system starting in the bigger region approach the smaller one, 
while no further convergence guarantees can be given. 

In Section ^ we explain how the destabilizing effect of a given quantizer can be measured. This 
destabilization measure, in conjunction with an arbitrary stabilizing feedback law and a corresponding 
Lyapunov function, can be used to determine an ultimate bound on solutions. One example of such a 
destabilization measure is the worst-case quantization error max x& x> \q(x) — x\. However, it turns out 
that there exist other destabilization measures which are actually more suitable in the present context. 
Although the parameters of the control system are used in the stability analysis, the destabilization 
measure itself is a function of the quantization regions and quantization points only. The quantizer 
design problem then naturally reduces to an optimization problem which consists in minimizing such a 
measure over all quantizers satisfying the information constraint. 

After casting quantizer design as an optimization problem, we proceed to explain how techniques 
from (optimal) facility location yield new insights into this problem as well as efficient algorithms for 
solving it. Facility location problems concern the location of a fixed number of facilities that provide 
service demanded by users; the objective is to minimize the average or maximal distance from sets of 
demand points to facilities. We focus here on settings continuous in the location of both the facilities 
and the demand points (both facilities and demand points take values in a continuum of points, e.g., in a 
polytope or ellipsoid). Facility location problems are surveyed in Q. Computational geometric aspects 
in continuous facility location are discussed in |l9|, |l^] and indirectly in textbooks on computational 
geometry Relevant background on computational geometric methods in locational optimization is 
provided in Section ^. 

For example, a classical problem of interest in locational optimization is the so-called multimedian 
problem. It consists in choosing a collection of N points qi, q2, ■ ■ ■ , qn in a bounded region V C W 1 
so as to minimize the quantity E(mm ie n ...m 1% — x \ 2 )i where the expected value is computed with 
respect to some probability density function on T> and | • | denotes the Euclidean norm. A solution of 
this problem is given by centroidal Voronoi tessellations] see [^, p|. Within the context of quantization 
and information theory, the multimedian problem is known as the classic fixed-rate minimum-distorsion 
quantizer design M. One of the early references on this problem is the work by Lloyd [13], who obtains 



optimality conditions and introduces a famous insightful algorithm. The multimedian problem is related 
to the problem of state moment stabilization of linear systems with limited data rate |l6| . 

Since we are working in the deterministic setting, we will find that the problem relevant for our 



purposes is the multicenter problem, discussed in [23, 22]. This is a somewhat less frequently encountered 



variant of the multimedian problem, which is obtained by replacing the expected value by the worst-case 
value; it can also be stated as the problem of covering a given region with overlapping balls of minimal 
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radius. The connection between the quantized control problem and the multicenter problem, although 
very natural, apparently has not been pursued before. In Section || we present a general formulation 
of the multicenter problem with weighting factors. We then discuss solutions of several versions of this 
problem which arise in the present context, all in terms of suitable Voronoi tessellations. We show 
how existing algorithms can handle some of the more apparent approaches, and then develop a novel 
algorithm which gives less conservative results. We note that while weighted multimedian problems are 
commonly encountered, our formulation is novel in that it introduces weighting factors in the multicenter 
problem. 

Simulation studies that we are currently conducting (as well as existing studies of the related multi- 
median problem, such as |gj) indicate that by solving the quantized feedback stabilization problem with 
the help of locational optimization techniques, one may obtain quite interesting quantization patterns. 
In particular, for planar systems a typical quantization region is a hexagon. Consequently, hexagonal 
quantization regions are capable of achieving better performance than more traditional rectangular ones. 

2 Quantization and stability 

We assume that the stabilization problem in the absence of quantization has been solved, in the sense 
that a state feedback control law is known such that the origin is a globally asymptotically stable equi- 
librium point of the ideal closed-loop system. In the presence of quantization, we adopt the "certainty 
equivalence" control paradigm; namely, we let the same control law act on the quantized state q(x). The 
problem under consideration is to characterize the destabilizing effect of a quantizer q, with the goal of 
obtaining an ultimate bound on solutions of the closed-loop system starting in a given bounded region. 
We first discuss this problem for general nonlinear systems and then develop more specific results for 
linear systems, moving from simpler but conservative to more complicated but sharper formulations. 

2.1 Nonlinear systems 

We start with the general situation where the process to be controlled is modeled by the nonlinear 
system 

x = f(x,u), xeR n , u £ R m . (1) 

All vector fields and control laws are understood to be sufficiently regular (e.g., C 1 ) so that existence and 
uniqueness of solutions are ensured. Suppose that some nominal static feedback law u = k(x) is given 
(with minor changes, dynamic feedback laws can also be used). In the presence of state quantization, 
we consider the feedback law u = k(q(x)) and the corresponding closed-loop system 

x = f(x,k(q(x)))=f(x,k(x + e)) (2) 

where 

e := q(x) — x 

represents the quantization error. 

Besides stabilizing the nominal system ([!]), the feedback law k clearly must possess some robustness 
property with respect to the measurement error e. To this end we impose the following assumption: 
there exists a C 1 function V : W 1 — > R such that for some class /Coo functions^ ai, a.2, CK3, p and for all 

2 Recall that a function a : [0, 00) — > [0, 00) is said to be of class K. if it is continuous, strictly increasing, and at(0) = 0. 
If a 6 K, is unbounded, then it is said to be of class /Coo. A function j3 : [0, 00) x [0, 00) — > [0, 00) is said to be of class K.C 
if f3(-, t) is of class /C for each fixed t > and /3(s, t) decreases to as t — > 00 for each fixed s > 0. 
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x, e G W 1 we have 

ai(M) < ^(x) < a 2 (M) 

and 

W > p(|e|) — /(x,fc(x + e)) < -a 3 (|s|). 

(Here and below, | • | denotes the standard Euclidean norm.) This amounts to the property that the 
control law u = k(x) input-to- state stabilizes the closed- loop system with respect to the measurement 



error e [20, The above assumption is rather restrictive and can be relaxed at the expense of 

obtaining weaker results (for linear systems, however, it is an automatic consequence of closed-loop 
asymptotic stability for e = 0). There is also considerable research on designing control laws satisfying 
this assumption. These issues are discussed elsewhere [|lT[ [llj . 

Pick a positive number M and consider the ball Bm '■= { x G M n : \x\ < M}. Consider the worst-case 
quantization error 

A := max |e| (3) 

xeB m 

(this quantity is sometimes also referred to as the sensitivity of the quantizer). The following result is 
then straightforward to prove (see (T^[). 



Lemma 1 Assume that 

Then the sets 

and 



ai(M) > a 2 op(A). (4) 
Hi := {x€R n : V{x) < ai(M)} (5) 
K 2 := V(x) <a 2 o p(A)} 



are invariant regions for the system (|2|). Moreover, all solutions of (|2|) that start in the set 1Z\ enter 
the smaller set 1Z 2 in finite time. An upper bound on this time is 

= ati(M) - a 2 o p(A) , . 

a 3 op(A) ' 1 ' 

This lemma implies, in particular, that all solutions starting in TZi at time t = to satisfy the ultimate 
bound 

\x(t)\ < a^ 1 oa 2 op(A) Vt>t + T (7) 

with T given by the formula (||). We regard the quantity A defined by (||) as a destabilization measure 
of the quantizer q. For given feedback law k and Lyapunov function V, an ultimate bound on solutions 
can be described by a class /Coo function of this measure as shown above. It is not hard to see that 
if the number N of quantization regions is sufficiently large, then A can be made small enough for 
the inequality (||) to hold. Minimizing A — and consequently the size of the attracting invariant region 
1Z\ — over all possible choices of the quantizer q corresponds to the following optimization problem: 

min max max |o~ — x\ (8) 
Q,W ie{i,...,N} xeWi 

where Q = {q±, . . . , q^} is a set of quantization points and W = {W±, . . . , Wn} is a partition of Bm 
into quantization regions. (We could work with partitions of TZ± rather than Bm, but this requires the 
knowledge of P and also may be less computationally feasible for non-quadratic V.) The optimization 
problem <M) is known as the multicenter problem in computational geometry; we defer its detailed 



discussion until Section 3.1. 
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2.2 Linear systems 

We now specialize to the case when the process is described by the linear system 

x = Ax + Bu, x G M. n , u G R m . (9) 

The linear system structure can be utilized to define a less conservative destabilization measure. Suppose 
that the system (|9|) is stabilizable, so that for some matrix K the eigenvalues of A + BK have negative 
real parts. Then there exists a unique positive definite symmetric matrix P such that 

(A + BK) T P + P(A + BK) = -I. (10) 

We let A m i n (-P) and A max (P) denote the smallest and the largest eigenvalue of P, respectively. We 
assume that B / and this is no loss of generality because the case of interest is when A is not 

a stable matrix. 

The quantized state feedback control law 

u = Kq(x) 

yields the closed-loop system 

x = Ax + BKq(x) = (A + BK)x + BKe (11) 
where e := q{x) — x is the quantization error as before. The derivative of the function 

V(x) := x T Px 

along solutions of the system (|Tl]) satisfies 

V = - x T x + 2x T PBKe < -\x\ 2 + 2\x\\PBKe\. (12) 
For an arbitrary small e > 0, we can rewrite this as 

V = ^M 2 - — - — -Id 2 + 2(1 + e)\PBKe\ 2 - I - 1 \x\ - J 2(1 + e)\PBKe\ ) 

l + e l 1 2(l+e)' 1 V n 1 \^2(lT^ 1 1 v v ;| 

< —\x\ 2 - — r f \x\ 2 - f 2(1 + e)\PBKe 

~ l + e l 1 2(1 + e) V V 

Therefore, we have 

\x\ > 2(1 + 6)1^5^61 =► ^<~y^I x I 2 - ( 13 ) 

Take a positive number M. There are two ways to proceed from here. One is to use the inequalities (O) 
and \PBKe\ < \\PBK\\\e\, where || • || stands for the induced matrix norm, to obtain 

|z| >2(l+e)||PB.K'|||e| => V < -^^l^l 2 - ( 14 ) 

Consider the ellipsoids 

Tlx := {x G R n : x T Px < A min (P)M 2 } (15) 

and 

1Z 2 := {x G R n : x T Px < A max (P)4(l + e) 2 \\PBK\\ 2 A 2 } (16) 

where A is the worst-case quantization error defined by (Q). Then we have the following linear coun- 
terpart of Lemma |l|. 



Lemma 2 Assume that 



A min (P)M z > A max (P)4(l +e)' z \\PBK\\ 2 A 2 . (17) 



Then the ellipsoids 1Z\ andlZ2 are invariant regions for the system (11). Moreover, all solutions of (p"l|) 
that start in the ellipsoid 1Z\ enter the smaller ellipsoid IZ2 in finite time. An upper bound on this time 
%s 

m A min (P)M 2 - A max (P)4(l + e) 2 \\PBK\\ 2 A 2 



A\\PBK\\ 2 A 2 (l + e)e 
Consequently, an ultimate bound on solutions starting in 1Z\ at time t = to is 



(18) 



\x(t)\<J^§-2(l+e)\\PBK\\A Vt > t + T 



with T given by the formula (jig) . Decreasing e to 0, we see that solutions (asymptotically) approach 
the ellipsoid 

{x G M n : x T Px < A max (P)4||PBiY || 2 A 2 }. 

Thus we still consider A as a destabilization measure. As in the nonlinear setting, this leads to the opti- 
mization problem (||). If N is large enough, then A can be made small enough so that the inequality ([HI) 
holds (for a given feedback gain K). 

Another approach is to work with (|D|) directly, avoiding the use of the induced norm ||PBif||. 
Define 

^■pbk '■= max \PBKe\. 

The result of Lemma |2| still holds if ||Pi?i^|| 2 A 2 is replaced by A 2 PBK everywhere in the statement of 
that lemma. This yields a less conservative ultimate bound and motivates the following optimization 
problem: 

min max max \PBK(qi — x)\ (19) 
Q,W ie{i,...,N} xeWi 

where Q is a set of quantization points as before and W is a partition of 1Z\ into N regions. This problem 
is in general lower-dimensional compared to (Jsj) because the subspace ker(PBK) can be ignored. (Note 
that PBK is a singular matrix whenever m < n.) Therefore, for the same N the optimal value for this 
problem will be significantly lower that for (|8|). However, A PBK is not really a destabilization measure 
in the sense used in this paper, because it depends on the feedback gain matrix K. While it gives better 
results for a fixed feedback law, quantizer design based on this destabilization measure needs to be 
redone if the feedback law is changed, and is not suitable for switching between several feedback laws. 
For these reasons we prefer to work with destabilization measures that are independent of a particular 
feedback law used. 

Remark 1 It is clear that in Lemma ||, the system's behavior is important only for x £ 1Z\ \ 7^2- We 
can thus work with partitions of 1Z\ \ IZ2 — or of a spherical annulus containing this set — rather than 
with partitions of 1Z\. This preserves or decreases the ultimate bound on solutions. □ 

Remark 2 Lemma |2| suggests that among stabilizing state feedback gains K, the ones that provide 
smaller ultimate bounds for the solutions of the quantized system are those with smaller values of the 
induced matrix norm ||PSif||, where P is given by (10). In this regard, it might be interesting to 



observe the following: if the open-loop system x = Ax is not asymptotically stable, then for every 
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stabilizing feedback gain K and the corresponding positive definite symmetric matrix P satisfying (ffO) 
we have 

\\PBK\\ > 1/2 (20) 
and the inequality is strict if x = Ax is unstable. To see this, use (|l2| ) and the definition of e to write 

V<-\x\ 2 (l-2\\PBK\\^-^ 

This formula will be used again several times in the sequel. Now, note that if q is chosen to take the 
value in a neighborhood of the origin, then the right-hand side of ( ^H) equals — |x| 2 (l — 2||Pfii^||) 
there, and so 1 — 2||PBif|| cannot be positive since A is not stable. It is also straightforward to show (p0|) 
directly: just multiply ( |10D on the left by v T and on the right by v, where v is a normalized eigenvector 
of A T P + PA with a nonnegative eigenvalue. □ 




2.3 Radial and spherical quantization 

In the above developments, the required bounds on the quantization error do not depend on the size 
of the state. This leads to uniform quantization, in the sense that quantization points are distributed 
uniformly over the region of interest. However, it is well known that more efficient quantization schemes 
are those which provide lower precision far away from the origin and higher precision close to the origin. 
Quantizers with a logarithmic scale are particularly useful; see J7j. Loosely speaking, with logarithmic 
quantization one has the same number of quantization points in the vicinity of every sphere centered 
at the origin in the state space, whereas with uniform quantization this number grows with the radius. 
This observation suggests introducing a "direct product" of one quantizer on a unit sphere and another 
along the radial direction, which is what we do next. 

Let us write 

x = |x|vers(x) 

where 

vers(x) := , — -. 

\x\ 

We represent the quantizer accordingly as 

q(x) = q r (\x\)q s (vers(x)) 

where q r takes N\ positive real values, q s takes N2 values on or inside the unit sphere, and N± and N2 
are some positive integers such that N1N2 < N. This means that we introduce two separate quantizers, 
one for \x\ and the other for vers(rc). Quantization points for the latter are allowed to be inside the unit 
sphere because, as will become clear later, requiring them to lie on the sphere would be too restrictive. 
(For example, for the case N2 = 1 it is easy to see that one quantization point at the origin achieves 
a smaller worst-case quantization error over the sphere than any point on the sphere.) The lattice of 
quantization points for the resulting overall quantizer q is formed by the N pairwise products of values 
of q r and q s . 

From the triangle inequality and the fact that |g s (vers(x))| < 1 for all x by construction, we obtain 

\q(x) — x\ < \q r '(\x\)q s ! (vers(x)) — |x|g s (vers(a;))| + | |x|g s (vers(x)) — |x|vers(x)| 
< — \x\\ + \x\ |g s (vers(x)) — vers(x)| 

+ I Q s (vers (x)) — vers (a;) | ) 



? r (M) 
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Using (|2l|) , we conclude that the derivative of V along solutions of the closed-loop system satisfies 



V < -Id 



1 - 2\\PBK\ 



+ k s (vers(x)) — vers(x) 



(22) 



Take two numbers A and e in the interval (0, 1). Then we have V < — e|x| 2 whenever the inequalities 



•(M) 



< 



A 



and 



are satisfied. 



|g s (vers(x)) — vers(x)| < 



2\\PBK\\ 

1- A-e 



2\\PBK\ 



(23) 
(24) 



Pick a positive number M. To handle (23), we take q r to be a logarithmic quantizer. Define 

A 



and 



a : = 1 



b:= 1 + 



2\\PBK\ 
A 



2\\PBK\\ 
noting from (|20|) that a > 0. Let 

g r (s) := (oV6 i-1 )Af for s G ((a/byM, (a/6) i_1 M) , i G {1, . . . , JVi} (25) 

and use arbitrary tie-breaking to define q r at the endpoints of the above intervals. Then it is easy to 
check that (|23|) holds for all x such that 

(a/b) Nl M < \x\ < M. 

In view of the condition (|24|), we introduce the worst-case quantization error on the unit sphere 
corresponding to q s : 

A s := max \q s (x) — x\. (26) 

\x\=l 

Performing a Lyapunov analysis similar to the one in the previous subsection, we obtain the following 
result. 

Lemma 3 Assume that the inequalities 

A min (P)Af 2 > A max (P)(a/6) Afl M 2 

and 



A, < 



1-A-e 
2\\PBK\\ 



(27) 



are satisfied. Then the ellipsoids 



111 ■= x T Px < A min (P)M 2 } 



and 



1Z 2 := {x£f: x T Px < X max (P)(a/b) Nl M 2 } 

are invariant regions for the system ( [TT|) . Moreover, all solutions of (|Tl| ) that start in the ellipsoid TZ\ 
enter the smaller ellipsoid IZ2 in finite time. An upper bound on this time is 

A min (P)M 2 - X max (P)(a/b) N ^M 2 
(a/b) 2N ^M 2 e 



For fixed N\ and N%, the quantity A s denned by ( fZq ) provides a destabilization measure (for q s ). 
When K is given and A s satisfies the inequality (p7|) for some A, e G (0, 1), we can construct q r via (p5|) 
and compute an ultimate bound on solutions using Lemma [| Minimizing A s corresponds to the 
following optimization problem: 

min max max \qf — x\ (28) 
2 s , w s ie{i,...,JV 2 } xew z s 

where Q s = {q{ , . . . , q% 2 } is a set of points on or inside the unit sphere and W s = {Wf, . . . , Wj^ 2 } is 
a partition of the unit sphere. An algorithm for solving this problem will be described in Section |3.1| . 
The quantity ([28|) will not exceed the right-hand side of (p7|) if N2 is sufficiently large. The values of 
iVi and N2 which yield the smallest ultimate bound seem to be difficult to compute analytically and in 
general depend on the stabilizing feedback gain K; however, for a given K we only have a finite number 
of choices for these integers and so can find optimal values by trying all of them. We remark that in the 
context of the multimedian problem, the idea of spherical coordinates quantization has been exploited 
before, and in particular the trade-off between the numbers of values for the radial and the spherical 



directions has been studied; see [24] and the references therein. 

It is straightforward to derive similar results using the norm defined by the Lyapunov function, 
i.e., ||x|| := V x T Px, instead of the Euclidean norm. This gives rise to an optimization problem on an 
ellipsoid rather than on a sphere. 



2.4 Radially weighted quantization 

The need for logarithmic quantization patterns is evidenced by the fact that it is the ratio |e|/|x|, and 
not the absolute value of the quantization error |e| itself, that needs to be small. This is clear from 
the formulas ( |l4|) and (^H) . The approach of Section leads to an "aligned" logarithmic quantization 
pattern, in the sense that quantization points on spheres of different radii are aligned along the same 
radial directions. This is a consequence of the fact that we treated radial and spherical quantization 
separately, having passed from the formula ( |2l| ) via the triangle inequality to the formula ( p2] ) and 
then to the two independent inequalities (|23| ) and (|24]). However, it is not hard to see that non-aligned 
quantization patterns can achieve better coverage. This suggests proceeding from ( |2l"| ) in a more direct 
fashion. 

To this end, pick two numbers M > m > and consider the ellipsoids 

Kx := {x£R n : x T Px < A min (P)M 2 } 

and 

TZ 2 := {i£R n : x T Px < A max (P)m 2 }. 

Define 

\q(x) — x\ , . 

A rw := max -. 29) 

x£7ei\7e 2 \x\ 



The following result then easily follows by virtue of (21). 
Lemma 4 Assume that 

A min (P)Af 2 > A max (P)m 2 (30) 

and 

A rw < „ ~ £ „ (31) 

~ 2 PM V ; 
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for some e > 0. Then the ellipsoids 1Z\ and IZ2 are invariant regions for the system (|ll|). Moreover, all 
solutions of (|ll|) that start in the ellipsoid H\ enter the smaller ellipsoid in finite time. An upper 
bound on this time is 

A min (P)M 2 - A max (P)m 2 



T 



m 2 e 



The quantity A rw defined by (^) provides another destabilization measure for q, in relation to a 
pair of numbers M > m > 0. Given a stabilizing feedback gain K, we can check the inequalities (p0|) 



and (31) and, if they are satisfied, obtain an ultimate bound on solutions from Lemma |4|. This leads us 



to the following optimization problem: 

min max max — : — - (32) 
Q,W ie{l,...,N} x&N l \x\ 

where Q = {q\, . . . , qjy} is a set of quantization points and W = {W\, . . . , Wn} is a partition of the 



annulus {x E M. n : m < \x\ < M} into quantization regions. The optimization problem (32) is different 
in structure from the ones we encountered earlier, and apparently has not been studied in the locational 
optimization literature. We henceforth call it the radially weighted multicenter problem. It turns out 
that while this problem is more challenging than the others, it is still computationally tractable. We 
will develop an algorithm for solving it in Section \5.2\ In particular, the inequality (^TJ) will hold for a 
given K if N is sufficiently large. 



2.5 Dynamic quantization 

So far we have been discussing quantizer design with the understanding that it can only be done once 
and cannot be changed on-line. The results obtained with this static quantization approach assert the 
existence of two nested invariant regions such that all trajectories of the quantized system starting in 
the bigger region approach the smaller one. Let us now suppose that we can recompute the locations 
of quantization regions and quantization points on-line. Then it is actually possible to achieve global 
asymptotic stability of the quantized closed-loop system. In this section we briefly explain how this 
can be done. The discussion closely parallels that of [||, [E3], and we refer the reader to these papers for 
details. In fact, the quantization schemes considered here all fit into the general framework developed 
in p|. 

A control strategy which incorporates dynamic quantization consists of two stages. The objective 
of the first, "zooming-out" stage is to find a number M such that x(to) belongs to the set 1Z\ defined 
by (|l5|) for linear systems or (||) for nonlinear systems. This is easy to accomplish, for example, by using 
a binary quantizer which takes one value inside the ball of a variable radius /i centered at the origin 
and another value outside this ball. If we set the control to zero and increase \x fast enough to dominate 
the expansiveness of the open-loop system, then at some time to we will know that the state lies inside 
the ball of radius /J.(to), at which point we can set M := y/ A max (-P) / A m in(-P) A^o) i n t ne linear case or 
M := a^ 1 o 02(^(^0)) i n the nonlinear case. 

During the second, "zooming-in" stage, we achieve asymptotic convergence to the origin by applying 
one of the results expressed by Lemmas |], ||, ||[ and |||, depending on which situation is being considered. 
We know that if an appropriate inequality involving the size of the quantization error is satisfied, then 
x(tQ + T) G 7^2- At time to + T, we redefine M so that IZi coincides with TZ\ for the new value of 
M. Recomputing the quantizer for this smaller region and repeating the procedure, we have x(t) — > 
as t — ► 00 (for linear systems, the convergence is exponential); one can also show that the origin is 
Lyapunov stable. This technique effectively mimics a logarithmic quantizer with a countable set of 
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values, of the kind considered in Q. It is important to recognize that since IZ2 is an invariant region 
for the closed-loop system, it is acceptable to use a conservative upper bound on T. 

When the quantization is static, our goal is to make the attracting invariant region IZ2 as small 
as possible. With dynamic quantization, we only need to ensure that IZ2 is a strict subset of 1Z\. 
(In fact, since convergence may become slower closer to the boundary of IZ2, we might want to make 
the annulus 1Z± \ IZ2 as small as possible, so that we can achieve better coverage on it with a given 
number of quantization points and obtain faster convergence.) The optimization problems formulated 
above remain relevant, except that by passing from static to dynamic quantization we basically pass 
from optimal to suboptimal quantizer design. In the context of Section |2.2| , this corresponds to simply 
making sure that A is small enough so that the inequality ( |17| ) is satisfied. In Section [2^ we just need to 
enforce the inequality (^7j) and can choose a small value for N\; for example, we can let N\ = 1, thereby 



eliminating radial quantization altogether. In Section 2.4 we do not need m to be small compared to 
M. 

Remark 3 In the above discussion we assumed that quantization regions and quantization points can 
be recomputed infinitely many times and as often as needed. Of course, having to recompute the 
quantization parameters too frequently may be undesirable. An interesting direction for future research 
is to incorporate a cost on the number and frequency of these "switching events" into the quantizer 
design problem. □ 



3 Continuous multicenter problems in facility location 

In this section we present a class of optimization problems related to the field of facility location; see 
the discussion in Section [l] and the survey ||. The facility location problem we consider will have as 
special cases the optimization problems studied in Section 0, and in particular the problems (||) from 



Section 2.1, (281) from Section 2.3, and (|32h from Section 



Let us review some preliminary concepts. Given a compact connected region T> C K n and a set of N 
points Q = {q\, . . . , qjy} in M n , the Voronoi partition V = {V±, . . . , V/v} of V generated by Q is defined 
according to 

Vi := {x £ V : \x — qi\ < \x — qj\ \/j 7^ i}. (33) 

When it is useful to emphasize the dependency on Q, we shall write V(Q) or Vi(Q). When V is a 
polytope in W 1 , each tyj is a polytope, otherwise Vi is the intersection between a polytope and V. The 
faces of the polytope which defines Vi are given by hyperplanes of points in W 1 that are equidistant 
from qi and qj, j 7^ i; among the latter, only "neighboring" points play a role. Note that this (standard) 
construction remains valid when T> is a lower-dimensional subset of M n , such as a sphere. We refer 
to H|, |l^] for comprehensive treatments of Voronoi partitions. 

Let Q = {qi, . . . , qpj} be a collection of points in M n and let W = {Wi, . . . , Wn} be a partition of 
T>. In what follows, we shall concern ourselves with the function 

H(Q,W):= max max <f>(x)f(\x - d) (34) 
ie{i,...,N} xeWi 

where eft : D — > [0, 00) is continuous non-negative and / : [0, 00) — ► [0, 00) is continuous, non-decreasing 
and unbounded. We also assume eft does not identically vanish on T>. We investigate the optimization 
problem 

minW(Q,W) (35) 
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and refer to it as the weighted multicenter problem. In general, TL is a nonlinear non-convex function of 
the locations Q and of the partition W. Accordingly, its global minima can be obtained only numerically 
via nonlinear programming algorithms. However, this and related facility location problems p3| , p^ , |6[ 
have some peculiar structure that helps us characterize optimal solutions and design useful iterative 
algorithms. Let us start by considering the weighted 1-center problem over T>, i.e., assume N = 1. 

Lemma 5 The function 7i\ : W 1 — > [0, oo) defined by 

Hi(q) := H({q}, {V}) = max<t>(x)f(\x - q\) 

is continuous, radially unbounded, and quasiconvex^ If f is convex and (f> is constant, then Tii is convex. 

Proof. The function 7i\ is continuous because it is the maximum of a compact family of continuous 
functions. Further, 7i\ is radially unbounded because for every x* G T> such that <f>(x*) > we 
have TL\{q) > <f)(x*)f(\x* — q\) and / is unbounded. To show the other statements, we invoke certain 
properties of convex and quasiconvex functions; see Sections 2.2 and 2.4 in Q. At fixed x, the function 
q i — ► f(\x — q\) is quasiconvex because it is the composition of a convex function with a non-decreasing 
function. Furthermore, if / is non-decreasing and convex, then q f(\x — q\) is convex because, at fixed 
x, it is the composition of a convex function with a convex non-decreasing function. If / is convex, then 
q i — ^ maxjgx) f(\x — q\) is convex because it is the pointwise supremum over a set of convex functions. 
For general 4> and /, the function H\ is quasiconvex because it is the weighted pointwise supremum of 
quasiconvex functions. □ 

Next, we let co(T>) denote the convex hull of T> and study global minima of TL\ as follows. 

Lemma 6 The set of global minimum points for 7i\ is compact, convex and has a non-empty intersec- 
tion with co(T>). If f is strictly increasing, then all global minimum points belong to co(T>). 

Proof. The fact that the set of global minimum points is compact and convex is an immediate 
consequence of continuity, radial unboundedness, and quasiconvexity. Let us prove the non-empty 
intersection with co(P). Suppose that q* G" co(D) is a global minimum point for H\. Let p* G co(T>) be 
the closest point to q* , i.e., p* := arguiin^^^ \q* — x\. Then \x — p*\ < \x — q*\ for all x G V, so that, 
for all x G V, we have 4>(x)f(\x — p*\) < 4>(x)f(\x — q*\) < max^gp (f)(x)f(\x — q*\) = Hi(q*). Therefore, 
Hi(p*) = rnax xG x> 4>( x )f(\ x ~ P*\) — an d P* a ^ so belongs to the set of global minimum points. 

When / is strictly increasing, the previous argument leads to Hi(p*) < 7~ii (q*), which contradicts the 
assumption that q* is a global minimum. □ 

The last two lemmas show that the weighted 1-center problem over T> is a quasiconvex optimization 
problem, i.e., it consists in minimizing the quasiconvex function Hi over the convex set co(T>). It 
is known that every quasiconvex optimization problem can be solved by iterative techniques (via a 
bisection algorithm solving a convex feasibility problem at each step; see Section 3.2 in Q). We call 
q*(D) a weighted center of the region D if it is a (possibly non- unique) global minimum point: 

q*(T>) := argmin max <j)(x)f(\x — q\). 
qeco(V) x&v 

Now, it is useful to consider the weighted multicenter problem again, and define W i— > Q* (W) as the 
map that associates to W a collection of iV (possibly non-unique) global minimum points for the cor- 
responding weighted 1-center problems; in other words, Q*{{W\, . . . , VFtv}) := {q*(W\), . . . , g*(M / 7v)}. 

3 Recall that a quasiconvex function is a function denned on a convex domain and with convex sublevel sets. 



12 



Note that these weighted centers are well defined in view of the above discussion since each W% is com- 
pact. Finally, define the Lloyd map (or the Lloyd algorithm) C : (Q, W) i— > (Q',W) where W' := V(Q) 
and Q 1 := Q*(W). 

Lemma 7 At fixed Q, the global minimum of W *—> TC(Q,W) is achieved at W = V(Q). At fixed W, 
a global minimum of Q i-> 7~i(Q, W) is achieved at Q = Q*(W). T/ie Lloyd map is a descent algorithm 
for the cost function H , i.e., an application of the map is guaranteed not to increase Tt. Given an initial 
pair (Q , Wo), the sequence {£ k (Qo, Wo), k > 0} approaches the largest set invariant under C such 
that H(C(Q,W)) = H(Q,W). 

Fixed points of the Lloyd map are pairs (Q,W) such that W is the Voronoi partition generated 
by Q and at the same time the points in Q are weighted centers for W. It is an open conjecture that 
global minima of TC are fixed points of the Lloyd map and that the iteration described in the lemma 
converges to such fixed points. Nevertheless, the algorithm is of interest to us because it is guaranteed 
to improve a given quantizer design and provides a good indication as to whether or not N is large 
enough to achieve the control objective. 

The classic Lloyd algorithm is tailored to the continuous multimedian problem as it appears, for 
example, in the problem of fixed-rate minimum-distorsion quantizer design; see |5[ 9]. The classic Lloyd 
algorithm differs from the one considered here only in the fact that the points in Q are moved to 
the centroids — as opposed to the weighted centers — of the respective Voronoi regions. (Centroids are 
solutions of the 1-median problems.) 

Next, we consider the specific settings that arise in the quantizer design problems discussed in the 
previous section. We study the multicenter problem (||), the spherical multicenter problem (pg|), and 
the radially weighted multicenter problem (|32|). To implement the Lloyd algorithm, two tasks must be 
carried out repeatedly. One consists in computing the Voronoi partition for a given set of points Q, 
which is accomplished by the standard procedure described earlier. The other amounts to computing a 
weighted center for each set Wi in a given partition. Thus for each of the specific multicenter problems 
studied below, we only need to explain how to solve the corresponding 1-center problem. Some additional 
remarks on the properties of these particular multicenter problems will also be provided. 



3.1 Multicenter problem 

Let us consider the problem (§) arising in Section |2.1| . The domain is a ball centered at the origin or, 
more generally, an ellipsoid, i.e., P = {i£ W 1 : x T Px < 1} for some positive definite symmetric matrix 
P. Note also that the problem ( |l9| ) arising in Section |2.2| reduces, via a linear change of coordinates, to 
the multicenter problem considered here in a lower dimension. 

In the problem (|8|), the weighting function (j) and the performance function / are the identity maps. 
Under these conditions, we refer to the optimization problem ( |35|) simply as the multicenter problem; 



see ]^J, 22]. The multicenter problem can be restated as the problem of covering the region T> with 
(possibly overlapping) balls of smallest radius. If B\ C M n is the unit ball centered at the origin, and if 
RB\ + q denotes the ball of radius R centered at a point q, the problem reads: 

mini? subject to U {RB\ + qA 5 L). 

ie{i,...,N} 

Let us analyze the 1-center problem. From Lemma || we know that this is a convex optimization 
problem. For each region Vi, the optimal solution q*(Vi) is the unique center of the minimal-radius 
enclosing sphere for Vi. When Vi C M 2 is a polygon, this sphere is referred to as the smallest enclosing 
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circle and algorithms are available to compute it; see Chapter 4]. When Vi C W 1 is a polytope, the 
smallest enclosing ellipsoid (in particular, sphere) can be computed via iterative convex optimization 
algorithms; see ji]. For a Voronoi region Vi near the boundary of T>, which is not a polytope, we can 
under-approximate it by a polytope generated by the vertices of Vi and suitable additional points on 
the intersection of Vi with the boundary of T>, and then compute the center of this polytope. For a 
sufficiently close under-approximation, this center will also be the center of Vi. 

When T> is a unit cube in M n , the optimal value of the problem (||) satisfies the bounds 



1 a <Jn 

< A < 



The upper bound is easily obtained by constructing a uniform cubical quantization pattern, while the 



lower bound is known as Sukharev's lower bound on dispersion [17, 15]. In the present case when T> is 
a ball, it is straightforward to obtain similar bounds by considering inscribed and superscribed cubes 
for T>. The upper bound can be used to evaluate the convergence of the Lloyd algorithm. When the 
lower bound on A is not small enough for the inequality (|j) or ( |i~7| ) to hold, it indicates that a different 
destabilization measure and/or a different stabilizing feedback law must be used. 

It is also useful to recall some known facts about the multimedian problem. It is conjectured in ||| 
that for TV sufficiently large, the optimal quantizer with respect to the uniform probability density is 
given by a tessellation (i.e., translation and rotation) of a fixed polytope, except near the boundary of 
the region of interest. In two dimensions, polygons that can give rise to such tessellations are equilateral 
triangles, rectangles, and regular hexagons. Among these, the hexagon is optimal, because it has the 
smallest mean-square quantization error with respect to its centroid per unit volume. This result remains 
true if we consider the worst-case rather than mean-square quantization error, which is the quantity 
being minimized in the multicenter problem. The hexagon achieves the smallest error with respect to 
its center. (For the unit volume regular hexagon this error is approximately 0.62, compared with 0.707 
for the square and 0.936 for the equilateral triangle; the unit-volume disk gives the error of 0.564 but 
disks cannot be used to obtain tessellations.) 



The spherical multicenter problem ( pq ) from Section corresponds to the setting where P = {i6 
W 1 : \x\ = 1} is the unit sphere in W 1 . Since the spherical multicenter problem is formulated in terms of 
the Euclidean distance in M n , Voronoi partitions of the sphere can be constructed as explained earlier 
for the general case. Voronoi regions will be spherical polytopes, i.e., intersections of polytopes with the 
unit sphere. The center of a spherical polytope Vi is the center of the minimal-radius enclosing sphere 
for Vi. We can consider a polytope in R n generated by the vertices of Vi and perhaps some other points 
in Vi. If enough points are taken, then the center of this polytope will also be the center of Vi. As we 
explained earlier, computing the center of a polytope is a computationally tractable task. 



3.2 Radially weighted multicenter problem 



Here, we study the problem (^) formulated in Section [2^, namely, 

\Qi - A 

mm max max — : — : — 
S,W ie{i,...,N} xeWi \x\ 

where the domain is the spherical annulus T> = {x € M. n : m < \x\ < M}. We consider the corresponding 
radially weighted 1-center problem over a set V GT>: 

\q-x\ . . 

mm max — : — : — . (So) 

q£co(V) XGV \x\ 
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The problem is well-posed because V is a subset of T> and therefore does not contain the origin. In 
what follows, we take V to be a polytope; if it is not, we approximate it by a polytope as before. We 
begin by making the following observation. 



Lemma 8 The optimal cost in the problem (|3q ) is smaller than 1 if and only if the set V is separated 
from the origin by a hyperplane. 

Proof. Suppose first that V is separated from the origin by a hyperplane, so that ^ co(V). Let q 
be the projection of the origin onto co(V), i.e., q := argmin,,. 6co /y) \x\. By construction, \x — q\ < \x\ 
for all x 6 V, hence max xg y \o~ ~ x \/\ x \ < 1- This implies that the optimal cost in the problem (|36|) 
is less than 1. To prove the converse, suppose on the contrary that G co(V). This means that the 
origin lies on the line segment between two points x\,X2 £ V. For the optimal cost to be less than 1, 
the optimal point q* must belong to the open ball {q S W 1 : \q — x\\ < \xi\} as well as to the open 
ball {q € M. n : \q — xg| < l x 2|}- But the intersection between these two sets is empty which is a 
contradiction. □ 

We shall henceforth assume that the set V is separated from the origin by a hyperplane. For N 
sufficiently large, the initial quantization points can be chosen in such a way that each of the resulting 
Voronoi regions indeed has this property. Since the Lloyd algorithm decreases the cost, Lemma |8] implies 
that all Voronoi regions will then have this property at every step of the iteration. 

From Lemma [| and ^ we know that the problem (^6|) is quasiconvex and can thus be handled by 
iterative convex optimization algorithms. However, its structure can be utilized to obtain a solution 
more constructively. Let us first present an equivalent formulation of this optimization problem. 

Lemma 9 Let V be a polytope separated from the origin by a hyperplane. Consider the problem of 
finding the sphere with center c and radius r which encloses V and minimizes r/\c\. Let (c*,r*) be the 
parameters of the optimal sphere. Then the optimal value for the problem (36) is 7* := r*/|c*| and the 
optimal point is q* : = (l — (7*) 2 )c*. 



Proof. Let f q (x) := \q — x\/\x\. In the problem (pq), we search for q that minimizes the value of the 
function f q on its smallest level set enclosing V. For any 7 > 0, the 7-level set of f q is described by 

I |2 2i |2 n 

\q — X\ — 7 \X\ = L). 

Because V is separated from the origin by a hyperplane, we know from Lemma || that the optimal value 
of 7 is smaller than 1. Thus from here on we will only be interested in 7 < 1. A square completion 
argument leads to 



1 



1-7 



2 



2| |2\ 
7 |x| ) 



1-7 



(1 - 7 2 ) 2 



so that the 7-level set of f q is the sphere \x — c| 2 = r 2 , with center c := qj (1 — 7 2 ) and radius r := r y\c\. 
In the new variables (c, r), we must minimize 7 = r/\c\ among all spheres enclosing V . □ 

Note that the point q* belongs to co(V) by Lemma ^, while c* might not. The above result leads 
us to considering the problem 

r 2 

min 7 2 (c, r) := --77 subject to |c — Vi\ 2 <r 2 , i = l,...,p (37) 

rel cr 



where Vi, . . . , v p are the vertices of the polytope V . This is an optimization problem subject to inequality 
constraints, which can be solved with a finite number of computations. The idea is to enumerate active 
constraints, according to the procedure described in the following algorithm: 
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1: for all subsets S of the set of vertices of V do 

2: compute the (eg, rs)-sphere minimizing 7 2 among all (c, r)-spheres touching all points in S 
3: end for 

4: discard all (cs, rs)-spheres not containing all vertices of V 

5: find global minima for ([}?]) by comparing the values of r$/\cs\ 2 among all remaining candidate 
spheres 

Steps 4 and 5 are straightforward comparison checks. Regarding step 1, it turns out we can restrict 
our search to sets S containing at least two vertices of V, by virtue of the following result. 

Lemma 10 The optimal sphere for the problem (|3^) touches at least two vertices of V , i.e., at least 
two constraints are active at the minimum. 



Proof. The proof is by contradiction. Suppose that the optimal sphere touches only one vertex. We 
denote this vertex by v and assume, performing an affine coordinate change, that it has coordinates 
(1,0,..., 0) T . Let c = (xi, X2, • • • , x n ) T . Then we are led to minimizing 

2,- - (xi-i) 2 + Er= 2 ^ 2 1 , 1-2X1 

7 (xi,x 2 ,...,x n ) = — — = i + _ — — . (38) 



Let us show that this function has no critical points besides the pole at the origin and the zero at v. 
We have 

9 7 2 _2(x?-£IL2X 2 --il (39) 



and 



^ (£r=i^ 2 ) 2 

<9 7 2 _ -2x^(1 - 2xi) 

^" (Er=i- 2 ) 2 ' 



% ^ 1. (40) 



In view of the formula (|40|), every critical point satisfies either Xi = 1/2 or x, = for all i 7^ 1. In 
the first case, the formula (|39|) implies that we must have EIL2 %i = — 1/2, and this equation has no 
solution. In the second case, ([^) gives two solutions: xi = (pole at 0) and x\ = 1 (zero at u). The 
pole at the origin is not a minimum. The zero at v corresponds to the sphere of radius centered at 
v, which is not a feasible solution because it does not enclose V. In summary, we have shown that the 
optimal sphere cannot touch only a single vertex of V. □ 

Regarding step 2, we need to minimize 7 2 over spheres passing through two or more vertices of V. 
Spheres passing through I generic points in M n are parameterized by n + 1 — I variables. A convenient 
parameterization is obtained by intersecting hyperplanes of points equidistant from pairs of points from 
a given set. Coordinates of the points on the intersection are given by affine functions of n + 1 — I 
free parameters. Note that the radius r of the sphere is uniquely determined by its center c and the 
vertices of V which lie on the sphere. It is straightforward to verify that the function 7 2 in ( |37| ) is a 
rational function whose numerator and denominator are quadratic inhomogeneous polynomials in these 
free parameters, and that critical points of 7 2 are solutions of n + 1 — I quadratic equations in the same 
number of unknowns. According to Bezout's theorem, this genericaly gives 2 n+1 ~' candidate optimal 
spheres (see Q). Step 2 is completed by choosing the one with the smallest radius. 

As an example of step 2, let us work out the planar case. When n = 2, the problem reduces to 
finding critical points of 7 for circles passing through / vertices of V, where I > 1 by Lemma ||. Since 
for / > 2 there is at most one circle passing through the corresponding vertices, we only need to explain 
how to solve this problem for / = 2. For convenience, let us consider an affine change of coordinates 
which places the two vertices at (1,0) T and (— 1,0) T and the origin at some point (xo,yo) T - Without 
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loss of generality, assume that yo > 0. The center of the circle is denoted by c = (x, y) T . We know that 
c must be equidistant from the two vertices, hence x = 0. Then we have 



7 2 



and so 

9 7 2 _ 2y(x 2 + (y- yo) 2 ) - 2(y - y )(l + f) _ -2y 2 y + 2y(x 2 + y 2 - 1) + 2y 



°y (x 2 + (y-y ) 2 ) 2 (x 2 + (y-y ) 2 ) 

Equating the numerator to 0, we arrive at the equation 

-y 2 yo + y(xl + y\ - l) + y = o. 

In the special case when yo = 0, this reduces to — 1) = 0. Since xq = ±1 corresponds to one of the 
vertices being at the origin, which cannot happen by our earlier assumption, the solution is y = (as 
is also clear from symmetry). When yo 0> the minimum is achieved at 



x 2 + y 2 - 1 - V(x 2 +y 2 - l) 2 +4y 2 
(Note that this goes to as yo approaches or oo.) 
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